{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0c3a5122",
   "metadata": {},
   "source": [
    "## Introduction to signal processing: FFT\n",
    "\n",
    "Highly recommended to watch the youtube video by 3Blue1Brown to understand the concepts of Fourier Transform.\n",
    "\n",
    "[![IMAGE ALT TEXT HERE](https://img.youtube.com/vi/spUNpyF58BY/0.jpg)](https://www.youtube.com/watch?v=spUNpyF58BY)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03c4333d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcaff8e1",
   "metadata": {},
   "source": [
    "## Toy Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc1b664d",
   "metadata": {},
   "source": [
    "#### Ex 1: Create a sinosoidal signal with frequency 1 Hz and sample it at 10 Hz. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b62bb047",
   "metadata": {},
   "outputs": [],
   "source": [
    "# simple sinus signal\n",
    "fs = 10 # sampling frequency\n",
    "t = ?\n",
    "s = ?\n",
    "plt.figure(figsize=(8, 2))\n",
    "plt.plot(?, ?)\n",
    "plt.xlabel('Time [s]')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2241446e",
   "metadata": {},
   "source": [
    "#### Ex2: Compute the FFT of the signal and plot the magnitude of the FFT.\n",
    "**Hint:**\n",
    "The output of the FFT is a complex number. The magnitude of the FFT is given by the absolute value of the complex number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "898dde45",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the frequency spectrum\n",
    "s_fft = \n",
    "\n",
    "freq = ?\n",
    "\n",
    "# plot the raw signal and output of the FFT\n",
    "fig, axes = plt.subplots(3, 1, figsize=(8, 6), sharex=True)\n",
    "axes[0].plot...\n",
    "\n",
    "plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b12fb34e",
   "metadata": {},
   "source": [
    "#### Ex 3: How does the frequency spectrum look like if the signal does not center around 0? Why?\n",
    "\n",
    "**Hint:**\n",
    "\n",
    "$$s(f) = \\int_{-\\infty}^{\\infty} s(t) e^{-2\\pi ift} dt$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dea28e7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# simple sinus signal\n",
    "fs = 10 # sampling frequency\n",
    "t = ?\n",
    "s = ?\n",
    "f, axes = plt.subplots(2, 1, figsize=(8, 4))\n",
    "axes[0].plot ... \n",
    "axes[0].set_xlabel('Time [s]')\n",
    "axes[0].set_ylabel('Amplitude')\n",
    "\n",
    "# plot the frequency spectrum\n",
    "s_fft = ?\n",
    "\n",
    "freq = ?\n",
    "axes[1].plot...\n",
    "axes[1].set_xlabel('Frequency [Hz]')\n",
    "axes[1].set_ylabel('Amplitude')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aba64e2d",
   "metadata": {},
   "source": [
    "### Ex 4: Can the real part of the FFT output be negative? Why? Modify the signal above to test your hypothesis.\n",
    "\n",
    "**Hints**:\n",
    "$$e^{i \\theta} = \\cos(\\theta) + i \\sin(\\theta)$$\n",
    "$$sin(2 \\pi f t) = (e^{2\\pi ift} - e^{-2\\pi ift}) / (2i)$$\n",
    "$$cos(2 \\pi f t) = (e^{2\\pi ift} + e^{-2\\pi ift}) / 2 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3002e800",
   "metadata": {},
   "source": [
    "Phase shift in time domain corresponds to multiplication by a complex exponential in the frequency domain\n",
    "$$s(f) \\xrightarrow{\\text{time shift by } \\theta} s(f)e^{-i2\\pi f\\theta}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cea14e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# simple sinus signal\n",
    "fs = 10 # sampling frequency\n",
    "t = ?\n",
    "s = ?\n",
    "f, axes = plt.subplots(2, 1, figsize=(8, 4))\n",
    "axes[0].plot...\n",
    "axes[0].set_xlabel('Time [s]')\n",
    "axes[0].set_ylabel('Amplitude')\n",
    "\n",
    "# plot the frequency spectrum\n",
    "s_fft = ...\n",
    "\n",
    "freq = ...\n",
    "axes[1].plot...\n",
    "axes[1].set_xlabel('Frequency [Hz]')\n",
    "axes[1].set_ylabel('Amplitude')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0c8774b",
   "metadata": {},
   "source": [
    "### Ex 5: Can you modify the signal so that its frequency spectrum only has non-zero values for positive frequencies?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b97352f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# simple sinus signal\n",
    "fs = 10 # sampling frequency\n",
    "t = ...\n",
    "s = ...\n",
    "\n",
    "f, axes = plt.subplots(2, 1, figsize=(8, 4))\n",
    "axes[0].plot...\n",
    "axes[0].set_xlabel('Time [s]')\n",
    "axes[0].set_ylabel('Amplitude')\n",
    "\n",
    "# plot the frequency spectrum\n",
    "s_fft = ...\n",
    "\n",
    "freq = ...\n",
    "axes[1].plot...\n",
    "axes[1].set_xlabel('Frequency [Hz]')\n",
    "axes[1].set_ylabel('Amplitude')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38533a22",
   "metadata": {},
   "source": [
    "## Real world acoustic signal from a railway track \n",
    "- sampled at 20 kHz\n",
    "- an acoustic signal "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97aadde7",
   "metadata": {},
   "outputs": [],
   "source": [
    "signal_file_path = 'acoustic_signal.npy'\n",
    "x = np.load(signal_file_path)\n",
    "Fs = 20000\n",
    "x.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf85408c",
   "metadata": {},
   "source": [
    "#### Ex 5: extract accoustic data sample from the mortar class under speed 80, plot the signal and its FFT."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74f54f0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = ... # compute the FFT\n",
    "Xfreq= ... # compute the frequency axis\n",
    "fig, axes = plt.subplots(2, 1, figsize=(8, 4))\n",
    "axes[0].plot(x)\n",
    "axes[0].set_xlabel('Time [s]')\n",
    "axes[1].plot(Xfreq[:len(X)//2], np.abs(X)[:len(X)//2])\n",
    "axes[1].set_xlabel('Frequency [Hz]')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5aba8e1d",
   "metadata": {},
   "source": [
    "#### Ex 6: To reduce noise, apply the provided `smooth` function and compare results with and without the smoothing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4cf1e11",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.signal\n",
    "\n",
    "def smooth(x, window_len=30, window='hanning'):\n",
    "    if x.ndim != 1:\n",
    "        x = x[:,0]\n",
    "    if x.size < window_len:\n",
    "        raise ValueError\n",
    "    if window_len<3:\n",
    "        return x\n",
    "    if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:\n",
    "        raise ValueError\n",
    "    s=x\n",
    "    #print(len(s))\n",
    "    if window == 'flat': # moving average\n",
    "        w=np.ones(window_len,'d')\n",
    "    else:\n",
    "        w=eval('np.'+window+'(window_len)')\n",
    "\n",
    "    #y=np.convolve(w/w.sum(),s,mode='valid')\n",
    "    y = scipy.signal.convolve(s,w/w.sum(),mode='same')\n",
    "    return y "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6a94ce7",
   "metadata": {},
   "outputs": [],
   "source": [
    "f, axes = plt.subplots(2, 1, figsize=(8, 4))\n",
    "X = ... # compute the FFT\n",
    "Xfreq= ... # compute the frequency axis\n",
    "axes[0].plot(Xfreq[:len(X)//2], X[:len(X)//2])\n",
    "axes[0].set_ylim(-100, 2000)\n",
    "axes[1].set_xlabel('Frequency [Hz]')\n",
    "\n",
    "X_smooth = smooth(X,500)\n",
    "Xfreq_smooth = ... # compute the frequency axis\n",
    "axes[1].plot(Xfreq_smooth[:len(X)//2], X_smooth[:len(X)//2])\n",
    "axes[1].set_ylim(-100, 750)\n",
    "axes[1].set_xlabel('Frequency [Hz]')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d027ff3",
   "metadata": {},
   "source": [
    "#### Ex. 7: We provide the code to plot 2 figures, one giving frequencies content of a signal from 0 to 10Khz and another from 0 to 1Khz. Why it is interresting to make zoom to the frequencies 0 to 1Khz?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7adfc88",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "X = np.abs(np.fft.fft(x))\n",
    "X = smooth(X, 500)\n",
    "Xfreq = ... # compute the frequency axis\n",
    "\n",
    "fig, axes = plt.subplots(2, 1, figsize=(8, 4))\n",
    "axes[0].plot(...)\n",
    "axes[0].set_xlabel(\"Frequency (Hz)\")\n",
    "\n",
    "axes[1].plot(...)\n",
    "axes[1].set_xlim([0,1000]) # zoom in into the first 1000 Hz\n",
    "axes[1].set_xlabel(\"Frequency (Hz)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adbfbe16",
   "metadata": {},
   "source": [
    "#### Ex 8: Subsample the signal x by N  Please explain why the spectrum amplitudes are decreasing when N increases, What can you do to counteract this effect?\n",
    "\n",
    "**Hint:**: use `scipy.signal.decimate()` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "199d5888",
   "metadata": {},
   "outputs": [],
   "source": [
    "# subsample the signal x by N\n",
    "N = ...\n",
    "x = np.load(signal_file_path)\n",
    "x_subsampled = ...\n",
    "x_subsampled.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f3f4bd6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the spectrum of the subsampled signal on top of the spectrum of the signal with max sampling frequency \n",
    "\n",
    "X = np.abs(np.fft.fft(x))\n",
    "X = smooth(X, 500)\n",
    "Xfreq = np.fft.fftfreq(len(X), 1/Fs)\n",
    "fig, axes = plt.subplots(2, 1, figsize=(8, 6))\n",
    "\n",
    "\n",
    "axes[0].plot(Xfreq[:len(Xfreq)//2], X[:len(Xfreq)//2], label=\"Original signal\")\n",
    "axes[0].set_xlabel(\"Frequency (Hz)\")\n",
    "axes[1].plot(Xfreq[:len(Xfreq)//2], X[:len(Xfreq)//2], label=\"Original signal\")\n",
    "axes[1].set_xlim([0,1000])\n",
    "axes[1].set_xlabel(\"Frequency (Hz)\")\n",
    "axes[0].legend()\n",
    "\n",
    "\n",
    "# Subsampled signal\n",
    "X_sub = np.abs(np.fft.fft(x_subsampled))\n",
    "X_sub = smooth(X_sub, 500)\n",
    "Xfreq_sub = np.fft.fftfreq(len(X_sub), 1/Fs)  # Modifications necessary here \n",
    "axes[0].plot(Xfreq_sub[:len(Xfreq_sub)//2], X_sub[:len(Xfreq_sub)//2], label='Subsampled signal') # Modifications necessary here \n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "axes[1].plot(Xfreq_sub[:len(Xfreq_sub)//2], X_sub[:len(Xfreq_sub)//2], label='Subsampled signal') # Modifications necessary here \n",
    "axes[1].set_xlim([0, 1000])\n",
    "axes[1].set_xlabel(\"Frequency (Hz)\")\n",
    "axes[1].legend()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62d91c61",
   "metadata": {},
   "source": [
    "#### Ex 9: What happen when N = 12 why it does not happen before when N=5?\n",
    "\n",
    "\n",
    "**Hint:**\n",
    "Take into consideration the new max frequency after subsampling considering the Shannon Nyquist theorem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "633cfa75",
   "metadata": {},
   "outputs": [],
   "source": [
    "# subsample the signal x by N\n",
    "N = 12\n",
    "x = np.load(signal_file_path)\n",
    "x_subsampled = ...\n",
    "\n",
    "X = np.abs(np.fft.fft(x))\n",
    "X = smooth(X, 500)\n",
    "Xfreq = np.fft.fftfreq(len(X), 1/Fs)\n",
    "fig, axes = plt.subplots(2, 1, figsize=(8, 6))\n",
    "\n",
    "\n",
    "axes[0].plot(Xfreq[:len(Xfreq)//2], X[:len(Xfreq)//2], label=\"Original signal\")\n",
    "axes[0].set_xticks(np.linspace(0,0.5*Fs,11),np.linspace(0,10000,11,dtype=int))\n",
    "axes[0].set_xlabel(\"Frequency (Hz)\")\n",
    "axes[1].plot(Xfreq[:len(Xfreq)//2], X[:len(Xfreq)//2], label=\"Original signal\")\n",
    "axes[1].set_xlim([0,1000])\n",
    "axes[1].set_xlabel(\"Frequency (Hz)\")\n",
    "axes[0].legend()\n",
    "\n",
    "\n",
    "# Subsampled signal\n",
    "X_sub = np.abs(np.fft.fft(x_subsampled))\n",
    "X_sub = smooth(X_sub, 500)\n",
    "Xfreq_sub = np.fft.fftfreq(len(X_sub), 1/Fs)\n",
    "axes[0].plot(Xfreq_sub[:len(Xfreq_sub)//2], X_sub[:len(Xfreq_sub)//2], label='Subsampled signal') # Modifications necessary here \n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "axes[1].plot(Xfreq_sub[:len(Xfreq_sub)//2], X_sub[:len(Xfreq_sub)//2], label='Subsampled signal') # Modifications necessary here \n",
    "axes[1].set_xlim([0, 1000])\n",
    "axes[1].set_xlabel(\"Frequency (Hz)\")\n",
    "axes[1].legend()\n",
    "plt.tight_layout()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}